############################################################################
# SE for the main time-invariant results
#
############################################################################

rm(list=ls())
graphics.off()

library(data.table)
library(fixest)
library(mvtnorm)
library(Hmisc)
library(xtable)

time0_se=Sys.time()

############################################################################
# RE results with SE
############################################################################

load("Data_new/RE_results.RData")
source("Code/0. CODE Auxiliary.R")
source("Code/12a. CODE RE for boot.R")

# Save point estimates
re_res = res_w2$table_mg[1:7]

# Bootstrap to obtain SE
dt[,ones:=1]
st = dt[,lapply(.SD,length),.SDcols="ones",by="state"]
state_list = as.numeric(st$state)
state_obs = as.numeric(st$ones)
ns = length(state_list)
B=1000
sseed=20250102

# Matrices to save bootstrap results
rboot = matrix(NA,nrow=7,ncol=1)
data_num = matrix(NA,nrow=4,ncol=B)
rownames(data_num) = c("states id now","states ori sampled","fips id now","fips ori sampled")

for (b in 1:B) {
  
  print(paste0("Sample number ",b))
  set.seed(sseed+b)
  
  ################################################
  ### Data
  ################################################
  
  state_now = sample(state_list,size=length(state_list),prob=state_obs,replace=TRUE)
  state_now = sort(state_now)
  dt_now = dt[state==state_now[1]] 
  dt_now[,state_new:=1]    
  for (k in 2:ns) {    
    dt_k = dt[state==state_now[k]] 
    dt_k[,state_new:=k]
    dt_now = rbind(dt_now,dt_k)
  }
  dt_now[,fips_new := state_new*100000+fips]
  data_num[1,b]=length(unique(dt_now$state_new))
  data_num[2,b]=length(unique(dt_now$state))
  data_num[3,b]=length(unique(dt_now$fips_new))
  data_num[4,b]=length(unique(dt_now$fips))
  dt_now[,state:=state_new]
  
  ################################################
  ### RE
  ################################################
  
  # FE estimates
  f = feols(yield ~ precr + ddr:factor(fips) | factor(fips) + factor(year)^factor(state) ,data=dt_now)
  s = summary(f,vcov="iid")
  dt_slopes= slopes_dt(s,names_var="_iid")
  setnames(dt_slopes,"coef_iid","coef")
  
  # Construct data at the county level
  dt_aux1 = dt_now[,lapply(.SD,sum),.SDcols="corn_area",by="fips"]
  setnames(dt_aux1,c("fips","w"))
  # Mean of temperature (weighted)
  dt_aux3 = dt_now[,lapply(.SD,wtd.mean,weights=corn_area),.SDcols="ddr",by="fips"]
  setnames(dt_aux3,c("fips","ddr_wmean"))
  # SD of temperature (weighted)
  dt_aux5 = dt_now[,lapply(.SD,wtd.var,weights=corn_area),.SDcols="ddr",by="fips"]
  dt_aux5[,ddr:=sqrt(ddr)]
  setnames(dt_aux5,c("fips","ddr_wsd"))
  # Merge 
  dim(dt_slopes)
  dt_slopes = merge(dt_slopes,dt_aux1,by="fips",all=TRUE)
  dt_slopes = merge(dt_slopes,dt_aux3,by="fips",all=TRUE)
  dt_slopes = merge(dt_slopes,dt_aux5,by="fips",all=TRUE)
  dim(dt_slopes)
  # Keep relevant variables and names
  dt_re=dt_slopes[,.(fips,coef,se_iid,ddr_wmean,ddr_wsd,w)]
  setnames(dt_re,c("se_iid","ddr_wmean","ddr_wsd"),c("se","z","sd_z"))
  setnames(dt_re,"fips","id")
  dt_re[,w0:=w]
  summary(dt_re)
  
  # Data at the county-year level with temperature and gamma to calculate marginal effects
  dt_ext = dt_now[,.SD,.SDcols=c("fips","year","corn_area","ddr")]
  setnames(dt_ext,c("id","t","w","ddr"))
  dim(dt_ext)
  summary(dt_ext)
  
  # RE
  res_w2 = RE_func_sim(data=dt_re,data_ext=dt_ext,min0=10,max0=200)

  # Add some summary variables
  rboot = cbind(rboot,res_w2$table_mg[1:7,])
  
  # Save 
  if (b%%100==0) {
    save.image("Data_new/Results_se_time_invariant.RData")
  }
  
}

### Checks on data_num
# Number of "new" states =31
summary(data_num[1,])
# Number of originally sampled states usually lower
summary(data_num[2,])
# Number of originally sampled fips usually lower than that "new fips"
table(data_num[3,]>data_num[4,])

# SE across bootstrapped samples
rboot = rboot[,-1]
re_se = apply(rboot,MARGIN=1,FUN=sd)

# Save
rm(list=setdiff(ls(),list("time0_se","re_res","re_se")))
list_keep = ls()
list_keep = c(list_keep,"list_keep","time1_se")

time1_se=Sys.time()
save.image("Data_new/Results_se_time_invariant.RData")

############################################################################
# Results with homogeneous coefficients and its se
############################################################################

rm(list=ls())
time0_se_p2 = Sys.time()
load("Data_new/Results_se_time_invariant.RData")
load("Data_new/Data_baseline.RData")

# Regression of yields on temperature and precipitation, with county FE and state-year FE
f = feols(yield ~  precr + ddr | factor(fips) + factor(year)^factor(state),data=dt, cluster="state")

hom_res = f$coeftable["ddr","Estimate"]
hom_se = f$coeftable["ddr","Std. Error"]

list_keep = c(list_keep,"hom_res","hom_se","dt","time0_se_p2")
rm(list=setdiff(ls(),list_keep))

############################################################################
# Bring GFE results and its se
############################################################################

load("Data_new/Data_for_CI.RData")

# Main estimates
gfe_res=wtd.quantile(dt_bounds$coef_gfe,probs=c(0.1,0.25,0.5,0.75,0.9),weights=dt_bounds$corn_area,na.rm=TRUE)
gfe_res=c(gfe_res,wtd.mean(dt_bounds$coef_gfe,weights=dt_bounds$corn_area,na.rm=TRUE))
gfe_res=c(gfe_res,wtd.var(dt_bounds$coef_gfe,weights=dt_bounds$corn_area,na.rm=TRUE))

# SE
gfe_se = apply(mg_gfe,1,sd)
gfe_se

list_keep = c(list_keep,"gfe_res","gfe_se")
rm(list=setdiff(ls(),list_keep))

############################################################################
# Bring FE and HPJ results
############################################################################

load("Data_new/Data_summary.RData")

# Percentiles, mean and variances: point estimates
data_all[,w:=corn_area]
data_all[is.na(w),w:=1]
name_vars = c("fe_het","mg_jk")
s=as.matrix(data_all[,lapply(.SD,wtd.quantile,weights=w,probs=c(0.1,0.25,0.5,0.75,0.9),na.rm=TRUE),
               .SDcols=name_vars])
s=rbind(s,as.matrix(data_all[,lapply(.SD,wtd.mean,weights=w,na.rm=TRUE),.SDcols=name_vars]))
s=rbind(s,as.matrix(data_all[,lapply(.SD,wtd.var,weights=w,na.rm=TRUE),.SDcols=name_vars]))
fe_res = s

#### Cross-sectional bootstrap of estimates of beta_i
dt = data_all[is.na(fe_het)==FALSE]
dt = dt[,.SD,.SDcols=c("state","fips","year","w","fe_het")]
summary(dt)
dt[,ones:=1]
st = dt[,lapply(.SD,length),.SDcols="ones",by="state"]
state_list = as.numeric(st$state)
state_obs = as.numeric(st$ones)
ns = length(state_list)
B=1000
sseed=20250102
name_vars="fe_het"

# Matrices to save bootstrap results
rboot = matrix(NA,nrow=7,ncol=1)

for (b in 1:B) {
  
  print(paste0("Sample number ",b))
  set.seed(sseed+b)
  
  ################################################
  ### Data
  ################################################
  
  state_now = sample(state_list,size=length(state_list),prob=state_obs,replace=TRUE)
  state_now = sort(state_now)
  dt_now = dt[state==state_now[1]] 
  dt_now[,state_new:=1]    
  for (k in 2:ns) {    
    dt_k = dt[state==state_now[k]] 
    dt_k[,state_new:=k]
    dt_now = rbind(dt_now,dt_k)
  }

  ################################################
  ### Results
  ################################################
  
  s=as.matrix(dt_now[,lapply(.SD,wtd.quantile,weights=w,probs=c(0.1,0.25,0.5,0.75,0.9),na.rm=TRUE),
                       .SDcols=name_vars])
  s=rbind(s,as.matrix(dt_now[,lapply(.SD,wtd.mean,weights=w,na.rm=TRUE),.SDcols=name_vars]))
  s=rbind(s,as.matrix(dt_now[,lapply(.SD,wtd.var,weights=w,na.rm=TRUE),.SDcols=name_vars]))
  rboot = cbind(rboot,s)
  
}

rboot = rboot[,-1]
fe_se = apply(rboot,MARGIN=1,FUN=sd)

list_keep = c(list_keep,"fe_res","fe_se")
rm(list=setdiff(ls(),list_keep))

time1_se_p2=Sys.time()
time1_se_p2-time0_se_p2

############################################################################
# Tabulate
############################################################################

hom = c(rbind(sprintf("%.3f",rep(hom_res,6)),paste0("(",round(hom_se,digits=3),")")),sprintf("%.3f",0),"")
fe = c(rbind(sprintf("%.3f",fe_res[,1]),paste0("(",round(fe_se,digits=3),")")))
jkk = c(rbind(sprintf("%.3f",fe_res[,2]),paste0("(",round(fe_se,digits=3),")")))
gfe = c(rbind(sprintf("%.3f",gfe_res),paste0("(",round(gfe_se,digits=3),")")))
re = c(rbind(sprintf("%.3f",re_res),paste0("(",round(re_se,digits=3),")")))

nn = c(paste0("Percentile ",c(10,25,50,75,90)),"Mean","Variance")
nnames = c(rbind(nn,rep("",length(gfe_res))))
df = data.table(nnames=nnames,hom=hom,fe=fe,jkk=jkk,gfe=gfe,re=re)
df

# Save in latex
l=dim(df)[2]
print(xtable(df,type="latex",display=c("s",rep("s",l))),hline.after = NULL,
      file="Results/Summary_table_time_inv.tex",include.rownames = FALSE,include.colnames=FALSE,
      sanitize.text.function = function(x){x},only.contents = TRUE)



